Bursting oscillations as well as the bifurcation mechanism in a non-smooth chaotic geomagnetic field model
Zhang Ran1, Peng Miao1, Zhang Zhengdi1, †, Bi Qinsheng2
Faculty of Science, Jiangsu University, Zhenjiang 212013, China
Faculty of Civil Engineering and Mechanics, Jiangsu University, Zhenjiang 212013, China

 

† Corresponding author. E-mail: 747524263@qq.com

Project supported by the National Natural Science Foundation of China (Grant No. 11472116), the Key Program of the National Natural Science Foundation of China (Grant No. 11632008), and the Postgraduate Research & Practice Innovation Program of Jiangsu Province, China (Grant No. KYCX17 1784).

Abstract

Based on the chaotic geomagnetic field model, a non-smooth factor is introduced to explore complex dynamical behaviors of a system with multiple time scales. By regarding the whole excitation term as a parameter, bifurcation sets are derived, which divide the generalized parameter space into several regions corresponding to different kinds of dynamic behaviors. Due to the existence of non-smooth factors, different types of bifurcations are presented in spiking states, such as grazing-sliding bifurcation and across-sliding bifurcation. In addition, the non-smooth fold bifurcation may lead to the appearance of a special quiescent state in the interface as well as a non-smooth homoclinic bifurcation phenomenon. Due to these bifurcation behaviors, a special transition between spiking and quiescent state can also occur.

1. Introduction

The theory of dynamical systems is an important tool for studying inherent mechanisms of various natural and human phenomena, especially when dealing with systems which contain complex nonlinear characters,[1] such as economic models,[2,3] nonlinear oscillator circuit models,[4,5] ecological models,[6] medical models,[7] and game theory.[8]

Recently, Gissinger proposed a model[9] which was used to study periodic or chaotic changes of sign of the amplitudes in dynamical systems, in addition to studying some specific situations of generator and hydrodynamic problems.[1] In the paper, he discussed the physical interpretation of this model and obtained several different dynamic regimes by choosing different values for the free parameters. Because the model offers a simple modelization of experimental and natural objects exhibiting complex dynamics, many researchers gave it further attention. For example, Elsonbaty et al. conducted bifurcation analysis of the model with one and two co-dimensions for earth magnetic field reversals.[11] Applying different transformations of the free parameters, the form is converted to another form to explore possible bifurcations, such as pitchfork bifurcation, Andronov–Hopf bifurcation[12,13] and homoclinic bifurcation,[14] which results in a critical value between the equilibrium points and bifurcation of this model. Many qualitative behaviors for certain lists of values of parameters are illustrated by bifurcation diagrams or LES spectra, demonstrating that some attractors can simultaneously coexist in the system’s behavior. Due to the special patterns of bursters in non-autonomous dynamical systems, accounting for the bursting oscillation mechanism[15] is an unavoidable problem for these types of flows with multiple time scales.

Based on previous research into the chaotic geomagnetic field model, a non-smooth factor and periodic external excitation are introduced to explore the effect of two time scales[16,17] in such Filippov systems.[18] In recent years, many researchers have explored piecewise smooth systems[19] with periodic excitation, which does not include a non-smooth factor. Two types of bifurcations will be investigated by exploring the relations between each group of parameters. In particular, the traditional method cannot be used to explore the complex dynamical behaviors on interfaces with unconventional bifurcation analysis. Because of the discontinuity of the vector field, the behaviors between two smooth subsystems and at their interface may be emphasized to show the effect of a non-smooth factor. A special transition from a quiescent state to a spiking state as well as a change in the attractor structure may occur when all parameters are fixed as a group at particular numbers. However, for non-autonomous vector fields such as a periodic excited oscillator, where there exists an order gap between the frequency of the excitation and the natural frequency, bursting phenomena can also be observed, the mechanisms of which still need to be investigated. In this letter, we focus on the effect of a non-smooth factor and periodic excitation for a chaotic geomagnetic field model.

2. Mathematical model

The chaotic[21] reversals system is a set of equations proposed as a model for reversals of a turbulent magnetic dynamo of a chaotic geomagnetic field. Gissinger introduced the new deterministic model as

where the parameters μ, ν, and γ are real. In fact, the amplitudes x and y can be regarded as the dipolar and quadrupolar components of the magnetic field. Noting that the amplitude z is coupling y and x in Eq. (1), we can therefore interpret z as a velocity mode which breaks the mirror symmetry of the earth core flow. The system (1) establishes the interaction between the dipolar and quadrupolar magnetic modes coupled by a mirror-antisymmetric velocity mode, and can be directly compared to geomagnetic data.

By applying different transformations to the system, it can be converted to a new form whose characteristic equation is presented in standard form. Then the center manifold theorem[21b] can be applied in order to reduce the dynamical system dimension and enable the study of the system around (0,0,0). The pitchfork bifurcation, Andronov–Hopf bifurcation and Horozov–Takens bifurcation can be observed by choosing suitable parameters and numerical simulation. Here, we introduce a periodic excitation and a non-smooth factor to model (2) and take suitable parameter values so that an order gap exists between the exciting frequency and the natural frequency. In order to explore the evolution of the dynamic behaviors of the new Filippov system with two time scales, when the geomagnetic field is disturbed by external factors, the mathematical model can be expressed as

where ω = A sin (Ωτ). When the exciting frequency is far less than the natural frequency, denoted by ω, i.e, ΩωN, the effect of two time scales may appear. Furthermore, the non-smooth boundary at x = 0, which divides the phase space into two regions corresponding to different subsystems, may affect the behavior of the dynamics.

3. Theoretical analysis
3.1. Equilibrium point analysis

Regular bifurcations may appear in the two smooth regions, while non-smooth bifurcation may take place on the non-smooth interface. When the exciting frequency is far smaller than the natural frequency, an order gap between the natural frequency and exciting frequency. Correspondingly, the system can be divided into two subsystems, expressed in the following forms: for x > 0,

and for x < 0,
For x = 0, only one equilibrium point can be observed, the stability of which is determined by the characteristic equation, expressed as
On the other hand, for x ≠ 0, the equilibrium point can be denoted by E(x,y,z) = (x,γx/vx2,γv/vx2). According to the Routh–Hurwitz criteria, a set of necessary and sufficient conditions for all roots of Eq. (5) to have a negative real part is given by Δ0 = a0 = 1 > 0, Δ2 = a1 = vμ + 1 > 0, and Δ2 = a1a2a3 > 0, where
and then the equilibrium point is asymptotically stable.

However, the number of equilibrium points of the autonomous system may change from one to three with variation of the parameters. The critical condition can be expressed as

where fold bifurcation may occur, corresponding to the jumping phenomenon between different equilibrium points. Similarly, Hopf bifurcation associated with a pair of purely imaginary eigenvalues may occur at
leading to possible periodic oscillations.

3.2. Unconventional bifurcation analysis

For system (2), an auxiliary parameter q ∈ [0,1] is introduced to discuss the behavior of the dynamics in the non-smooth interface Σ = {(x,y,z)|x = 0} through differential inclusion theory. The system (2) can be changed into

The dynamic evolution of the system in the region D (or D+) for q = 0 (or q = 1) is decided by the smooth sub-system F(X) (or F+(X)). However, the behavior which is analyzed according to different contact modes with the interface Σ = {(x,y,z)|x = 0} is complex.

The equation dx/dt = −yz + (2q − 1)a0 + w can be used to describe projection in the interface when the periodic excitation term exists. The interface

is defined to ascertain the relationship between scale and the non-smooth factor, where q is a continuous closed set. If a0 > 0, the supremum and infimum exist for any time τ = τ′, which can be represented respectively by
Therefore, two boundaries of the region at τR can be given by The corresponding equation Fs can be expressed as

When the parameter Ω is taken to be 0 < Ω = 0.005 ≪ 1, there are two obvious scales in the frequency domain between the natural frequency and the exciting frequency. In the next section, we focus on the mechanism of bursting oscillation with fixed set of parameters A = 5, Ω = 0.005, μ = 0.5, v = 3.5, γ = 0.1, and a0 = −1.

4. The mechanism of bursting oscillation
4.1. Bursting phenomenon

The phase portrait on the (x,y) plane as well as the time history related to the state variable x for v = 3.5 are plotted in Fig. 1. The projection is symmetric about the origin in the plane (x,y), moving in a counter-clockwise direction with period 2π/Ω. The trajectory alternates over time between spiking states and quiescent states to show periodic oscillations, which can be observed in the time history.

Fig. 1. (color online) Periodic bursting oscillation. (a) The phase diagram in (x,y). (b) The time process of the variable x. (c) Local magnification of the phase diagram near the origin. (d) and (e) Local magnification in panel (b). (f) Bursting oscillation mode in a periodic trajectory.

The vector field of system (2) is divided into two smooth areas by the non-smooth interface Σ: = {(x,y,z)|x = 0}. In Fig. 1(c), the trajectory contacts the interface several times during one periodic oscillation and shows obvious non-smooth characteristics. The connection between the trajectory and the non-smooth interface is mainly concentrated on three positions, represented by Ps0, Ps1, Ps2. In the region D, the trajectory enters the interface from the point Ps1 along the direction of the blue arrow, and continues moving until it arrives at the point Ps0 where the trajectory starts to enter the area D+. Correspondingly, the trajectory may move on the interface from the point Ps2 to the origin Ps0 with the direction of the red arrow in the area D+. In order to show clearly the non-smooth characteristics of the other two positions with violent oscillation, the corresponding position is amplified in Figs. 1(d)1(e). At the beginning of the oscillation, the trajectory passes through the interface between two regions. Then the trajectory touches but does not pass the interface with the amplitude reducing. When the amplitude is decreased a little, the trajectory no longer contacts the interface and only oscillates in the corresponding area.

It is necessary to point out the time when the trajectory is near the origin for much longer than other points of contact with the interface in Fig. 1(b). Therefore, a special quiescent state exists in interface around the origin for a periodic oscillation. In a periodic bursting oscillation, there is a special quiescent state in the interface near the origin. The oscillation presents four types of quiescent state and spiking state which are represented by QS±i (i = 1,2) and SP±i (i = 1,2), respectively, in Fig. 1(f).

4.2. Oscillation mechanism

Because two time scales and a non-smooth factor evolve in the vector field, a slow–fast effect, which typically appears in bursting oscillations in different forms, can be observed in the system. In order to explore the mechanism of these special dynamical behaviors, the transformed phase portrait is introduced.

As shown in Fig. 2, the trajectory starts to oscillate from the point A1, at which the slow parameter w takes the minimum value w = −5. With increasing time t, the trajectory may move almost strictly along the stable equilibrium branch EB−1 towards the right, which appears in the quiescent state QS−1 until it reaches the super-critical bifurcation point, producing a stable limit cycle LCS2 with unstable equilibrium branch EB−2. However, due to the slow passage effect, the trajectory may continue moving towards the right along the branch EB−2 for a period of time and still be in a quiescent state. The trajectory may enter into the spiking state SP−2 until it begins to gradually oscillate again around the branch EB−2. In the process of oscillation, the limit cycle LCS2 may disappear at point , where a non-smooth homoclinic bifurcation appears. At the same time, the saddle point E−0 (x,y,z) = (0,0,0.1), which belongs to the smooth subsystem F and sliding boundary ∂Σ, may also transform into the corresponding nodes Es0(y,z) = (0,0.1) of the sliding area Σs through non-smooth bifurcation.

Fig. 2. (color online) The bursting oscillation for v = 3.5. (a) The transition diagram in the plane (w,x). (b) The superposition of the transformation phase diagram and Fig. 5(b).

Before the trajectory reaches the point FB (w,x) = (−0.263, −1.606), it contacts the non-smooth interface at the point A2 (w,x) = (−0.305,0) which corresponds to the point (x,y,z) = (0,−0.01,0.162) of the traditional phase space. At this time, the corresponding auxiliary parameter q satisfies

Obviously, the auxiliary parameter q ∈ (0,1) is satisfied. Therefore, the point A2 is located in the sliding area Σs and the trajectory may remain on the non-smooth interface for a moment since it is governed by the planar system Fs. At the same time, the corresponding node Es0 (y,z) = (0,0.1) of the sliding area Σs will appear through the non-smooth fold bifurcation. Based on the form of Fs, an analytic expression for the projection is given by solving the ordinary differential equation
where y0 and z0 represent variable values at the initial time when the trajectory enters the sliding region Σs. The sliding boundary ∂Σ± may continue to be far away from the origin with increasing parameter w, making the equilibrium point Es0 pass the sliding boundary ∂Σ±. Therefore, the time at which the trajectory enters the interface from the point A2 and passes the sliding boundary ∂Σ± can be calculated by the values of the corresponding parameter w and is
It can be judged that there is sufficient time for the trajectory to converge to the smallest area of the node Es0 in the sliding region. This process can be called the sticky sliding state, which is considered to be a special quiescent state QSS1 on the non-smooth interface.

Fig. 5. (color online) The homoclinic orbit of the saddle E+0 at w = 1. (a) Local flow patterns near the homoclinic orbit. (b) The equilibrium branch with change of the slow variable and the related single parameter bifurcation of the generalized autonomous system.

Two important unconventional bifurcations may appear for different reasons with increasing values of the slow parameter w. One reason is that the node Es0 of the sliding boundary turns into the saddle E+0 of the smooth subsystem by the non-smooth fold bifurcation. Another reason is that the limit cycle LCS1 collides with the saddle point E+0 of the sliding boundary ∂Σ± by the non-smooth fold bifurcation, which can be briefly described as non-smooth homoclinic bifurcation.

After the parameter w crosses w = 1, the trajectory may oscillate and converge to the limit cycle LCS1 with homoclinic characteristic, leaving the non-smooth interface from the sliding boundary ∂Σ±. As shown in Fig. 3(a), relative to the non-smooth interface, the trajectory may exhibit an across-sliding oscillation mode in the interval w ∈ [1,2.061]. Since the sliding boundary continues to be further away from the origin, the characteristics of the limit cycle LCS1 are gradually weakened, and the oscillations of the trajectory appear to gradually intensify. The trajectory may exhibit across-sliding bifurcation at the point w = 2.061 and no longer pass the non-smooth interface, sliding on the interface until it reaches the point , shown in Fig. 3(b). At this moment, the limit cycle LCS1 and the non-smooth boundary ∂Σs+ are tangent to a point on the sliding boundary, where grazing-sliding bifurcation occurs. Then sliding boundary continues to stay away from the origin, and the limit cycle is no longer in contact with the non-smooth interface. Under the influence of the above behaviors, the trajectory leaves the sliding mode to follow the further evolution of the limit cycles LCS1 and be entirely in the smooth region D+ since it is driven by the smooth subsystem F+.

Fig. 3. (color online) Local large map of the spiking state SP+1. (a) w ∈ [1,2.061] and (b) w ∈ (2.061,2.605].

With further increase in the slow parameter w, the limit cycle LCS1 continues to converge and oscillate towards the unstable equilibrium branch EB+2 until it arrives at the super-critical Hopf bifurcation point HB+ (w,x) = (3.531,1.835), causing decreasing oscillation amplitude. The limit cycle may shrink to the equilibrium branch EB+2 and disappear after crossing the point HB+. Then the trajectory moves strictly along the stable equilibrium branch EB+1 until it arrives at the maximum point A4 (ω = 5). Due to the symmetry of the system, the oscillation mode will show the same dynamical behavior in another smooth region.

4.3. Analysis of three types of bifurcation

For slow variable w = 2.605, corresponding to in Fig. 2, the limit cycle LCS1 is tangent to the non-smooth interface and its minimum point (x,y,z) in the x-axial direction satisfies

Therefore the point P where the limit cycle is tangent to the non-smooth interface is exactly located at the sliding boundary ∂Σs+, shown in Fig. 2(b). With decreasing parameter w, the sliding boundary will move towards the origin, which indicates that the limit cycle will enter the sliding area Σs = {(y,z)|x = 0} when it returns to the interface with a counter-clockwise direction again. This solution corresponds to the grazing-sliding bifurcation of the generalized autonomous system (2) for w = 2.605.

In fact, the trajectory may slide on the interface for a period of time after detaching from the smooth area D+, and then return to the smooth area D+ to continue oscillating. The minimum value of its amplitude in the x-axial direction remains zero until another critical case, w = 2.061. The limit cycle LCS1 which corresponds to in Fig. 2 may slide in the interface Σ from point P1 (x1,y1,z1) to point P2 (x2,y2,z2), shown in Fig. 4(c) by the blue solid line, and return to the smooth area D+ at point P2. Through calculation, P1 and P2 satisfy

which means that the points P1 and P2 may be located at the sliding boundaries ∂Σs− and ∂Σs+, respectively. Once the slow parameter passes the point w = 2.061, the limit cycle LCS1 may reach the interface above the sliding boundary ∂Σs− and settle in the area Σ while the sliding boundary ∂Σ continues moving towards the origin, and then the cycle may enter the smooth area D directly through the interface.

Fig. 4. (color online) The stable limit cycle LCS1 with different values of the slow variable. (a) w = 2.8. (b) The critical case of edge-sliding at w = 2.605. (c) The critical case of cross-sliding at w = 2.061. (d) The mode of cross-sliding at w = 1.6.

The structure of the limit cycle is given in Fig. 5(b) for w = 1.6. The cycle may contact the interface again when it arrives at point P2, and then slides to the point P3 on the boundary ∂Σs+. Finally, the cycle may leave the interface and return to the smooth area. The oscillation structure of the limit cycle is obviously different from the case which is presented in the interval w ∈ (2.061,2.605). Therefore, it corresponds to the unconventional bifurcation, the across-sliding bifurcation for w = 2.061.

In particular, the homoclinic orbits[22] are produced around the unstable focus at the saddle point E+0 which belongs to the smooth area F+ and sliding boundary ∂Σs+ for w = 1, shown in Fig. 5(a), which causes the phenomenon that the cycle LCS1 may be cut off and disappear at the point E+0. Therefore, for w = 1, a special bifurcation called non-smooth homoclinic bifurcation[16] may happen when the cycle disappears at the point E+0. With decreasing parameter w, the saddle point E+0 may be changed into the node Es0 of the planar system by non-smooth fold bifurcation and the sliding boundary ∂Σ, including the saddle point E+0, may enter the sliding area Σs. In Fig. 5(b), the homoclinic orbit slide into the interface exhibiting a sliding-viscous phenomenon. On reaching parameter w = −1, the node Es0 may be transformed into the saddle E−0 of the smooth subsystem F.

5. Conclusion

The mechanism of bursting oscillation in a chaotic geomagnetic field model is described in detail by presenting relative diagrams and numerical calculation. It is found that conventional bifurcation can affect the attractor structure of bursting oscillations and may lead to transitions between spiking states and quiescent states, such as supercritical Hopf bifurcation. Besides that, the effect of non-smooth bifurcation cannot be ignored. As far as this article is concerned, we describe two aspects of the effect. Firstly, the grazing-sliding bifurcation and across-sliding bifurcation of the limit cycle LCS1,2 may bring two types of modes of oscillation, which are presented as the sliding mode and the sliding with crossing mode in spiking states. Thus the structure of the attractor is affected in bursting oscillation. Secondly, the non-smooth fold bifurcation appearing at the equilibrium point located at the sliding boundary ∂Σ± causes the special quiescent state QSs1,2 in the non-smooth interface. Meanwhile, the transformation between the spiking state SP±1 and the quiescent state QSs1,2 will be affected by the non-smooth homoclinic bifurcation which relates to fold bifurcation. Therefore, it is not sufficient to consider only the spiking state and quiescent state for bursting oscillation. The oscillation modes and possible bifurcations should be also analyzed in detail when the trajectory is located in the spiking state.

Reference
[1] Lei Q Jiang D 2014 J. Chin. Univ. Pet. 9 1415
[2] Liao X Li C Zhou S 2005 Chaos Soliton. Fract. 1 91
[3] Cai G Huang J 2007 Int. J. Nonlinear Sci. 3 213
[4] El-Sayed A M A Nour H M Elsaid A 2016 Appl. Math. Model. 5 3516
[5] Elsonbaty A R El-Sayed A M A 2017 Nonlinear Dynam. 2 1169
[6] Agiza H N Elabbasy E M El-Metwally H 2009 Nonlinear Anal. -Real 1 116
[7] Wang Q Y Lu Q S 2008 International Journal of Bifurcation & Chaos 4 1189
[8] Elsadany A A Agiza H N Elabbasy E M 2013 Appl. Math. Comput. 24 11110
[9] Gissinger C 2012 Eur. Phys. J.B 4 137
[10] Hughes D W Proctor M R E 1990 Nonlinearity 1 127
[11] Elsonbaty A Elsadany A A 2017 Chaos Soliton. Fract. 103 325
[12] Peng M Zhang Z Wang X 2017 Adv. Differ. Equ. 1 387
[13] Wang B Jian J 2010 Commun. Nonlinear Sci. 15 189
[14] Yu W Q Chen F Q 2010 Commun. Nonlinear Sci. 15 4007
[15] Sokoloff D D 2017 Izv. Phys. Solid Earth 6 855
[16] Li X H Bi Q S 2012 Chin. Phys. 21 060505
[17] Bi Q S Zhang Z D 2011 Phys. Lett. A 8 1183
[18] Morel C Morel J Y Danca M F 2017 Math. Comput. Simulat. 10 1
[19] Li S Liu C 2015 J. Math. Anal. Appl. 2 1354
[20] Han Q S Chen D Y Zhang H 2017 Chin. Phys. 26 128202
[21] Naschie M S E 2006 Chaos Soliton. Fract. 4 816
[22] Marciniak-Czochra A Mikelic A 2016 Vietnam Journal of Mathematics 1 1
[23] Tian R L Yang X J Cao Q J Wu Q L 2012 Chin. Phys. 21 020503